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ABSTRACT 

This paper studies the formation and evolution of binary supermassive black holes (SMBHs) in 
rotating galactic nuclei, focusing on the role of stellar dynamics. We present the first iV-body simu- 
lations that follow the evolution of the SMBHs from kiloparsec separations all the way to their final 
relativistic coalescence, and that can robustly be scaled to real galaxies. The A^-body code includes 
post-Newtonian (VM) corrections to the binary equations of motion up to order 2.5; we show that the 
evolution of the massive binary is only correctly reproduced if the conservative ITWand 27W terms 
are included. The orbital eccentricities of the massive binaries in our simulations are often found to 
remain large until shortly before coalescence. This directly affects not only their orbital evolution 
rates, but has important consequences as well for the gravitational waveforms emitted during the 
relativistic inspiral. We estimate gravitational wave amplitudes when the frequencies fall inside the 
band of the (planned) Laser Interferometer Space Antennae (LISA). We find significant contributions 
— well above the LISA sensitivity curve — from the higher-order harmonics. 

Subject headings: black hole physics - gravitational waves - galaxies: evolution - galaxies: interactions 
- galaxies: kinematics and dynamics - galaxies: nuclei 



1. INTRODUCTION 

Supermassive black holes (SMBHs) are commonly ob- 
served at the centers of nearby galaxies, and the existence 
of quasars at rcdshifts z w 6 implies that many of these 
SMBHs reached nearly their current masses at very early 
times. Bright galaxies are believed to form via hierarchi- 
cal merging of smaller galaxies, and galaxy mergers lead 
inevitably to the formation of binary SMBHs (Begelman, 
Blandford & Rccs 1980). At their formation, such SMBH 
binaries typically have separations a « Oh, where 

,,^ _^3.3pc^^^(^ -^^^j , (1) 

(see, e.g.. lMerritt fc Milosavlievicll2005f) . This "hard bi- 
nary separation" ah is defined in the standard way as 
the value of the binary semi-major axis at which passing 
stars arc ejected with velocities high enough to unbind 
them from the SMBHs; is a function of the binary re- 
duced mass fi = mi 771,2/ ("^i + "^2) and the ambient stel- 
lar velocity dispersion a {q = m2/mi < 1). At least one 
binary SMBH has probably been obser ved, with a pro - 
jected separation of about 7pc (Rodriguez et al.ll2006D . 
The quasi-periodic outbursts of the quasar OJ287 have 
also plausibly been modeled as a binary system with sep - 
aration ~ 0.05 pc and eccentricity ~ 0.6 ()Valtonenll2007D . 
A few other galactic systems are known to contain dual 
SMBHs with separations of a few k iloparsec, presumably 
the precursors of b inary SMBHs (jKomossa et all 120031 : 
iBianchi et a"ni2008[ ). 

Binary SMBHs may eventually coalesce, but only af- 
ter stellar- and/or gas-dynamical processes have first 
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brought the two SMBHs to separations small enough 
(~ 10"'^ pc) that gravitational wave emission is effec- 
tive. Whether Nature typically succeeds in overcoming 
this "final parsec problem," or whether long-lived binary 
SMBHs are the norm, is currently unknown. Persistence 
of the binaries would have a number of potentially im- 
portant consequences, including ejection of SMBHs from 
galaxy nuclei via three-body interactions and a reduc- 
tion in the mean ratio of SMBH m ass to galaxy mass 
(Volonteri, Haardt fc Madau 2003; iMadau et all 120041 : 
iLibeskind et al.|[2006f ). In addition, plans to detect grav- 
itational waves from the final inspiral of binary SMBHs 
by space-based observatories such the Laser Interferome- 
ter Space Antennae (LISA) (Sesana, Volonteri fc Haardt 
2007) might need to be re-considered if evolution of the 
binaries typically stalls at large separations. 

In a spherical, gas-free galaxy, the evolution of a bi- 
nary SMBH slows down at separations of ~ Oh because 
stars on orbits that intersect the binary - so-called "loss- 
cone orbits" - are depleted in just a few galaxy cross- 
ing times (e.g., Begelman et al. 1980, Makino et al. 
1993, Makino 1998, Berczik, Spurzem fc Merritt 2005; 
Merritt 2006). A massive binary can continue to evolve 
in such a galaxy only if the loss-cone orbits are some- 
how repopulated. Gravitational scattering is one possible 
mechanism for loss-cone repopulation, but it is only ef- 
fective in low-luminosity galaxie s with cent ral relaxation 
times below ~ 10 Gyr (jYu 2002; Merritt. Mi kkola fc- Szcll 
|2007() . Dense concentrations of gas can substantially ac- 
celerate the evolution of a massive b inary by increasing 
the drag on the indivi dual SMBHs (jEscala et all 12004 
120051: iDotti et al.ll2007ft . The plausibility of such gas ac- 
cumulations, with masses comparable to the masses of 
the SMBHs, is unclear however, particularly in the most 
massive galaxies. Galaxy merger simulations including 
gas have followed the two SMBHs until separations of a 
few tens of parsecs, rou ghly the hard bina ry separation 
(jKazantzidis ct al. 2005: iMaver et al.|[2007[ ). These sim- 
ulations contain useful information about the formation 
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and early evolution of SMBH binaries but - at least so 
far - do not have the resolution needed to address the 
final parsec problem. 

An alternative pathway exists for binary SMBHs to 
evolve beyond a « Oh, even in gas-free galaxies with 
long relaxation times. If the galaxy potential is non- 
axisymmetric, many of the stars will be on centrophilic 
orbits, i.e., orbits that pas s near the center of the g alaxy 
once per orbital period ([Gerhard fc Binn"evl Il985f ) . If 
even a small fraction of a galaxy's mass is associated 
with such orbits, the "feeding" rate of a central binary 
can be enormously enhanced compared w i th the rates in 
spheri cal galaxies (jMerritt fc Poonll2004l ). iBerczik et al.l 
(|2006[ ) explored this pathway in iV-body simulations con- 
taining particle numbers up to one million; the initial 
galaxy models were rotating, and in some cases the rota- 
tion was rapid enough to induce the formation of a triax- 
ial bar. Evolution of the SMBH binary was observed not 
to stall at a « ah in the triaxial models; furthermore the 
evolution rates exhibited no discernible dependence on 
particle number, as predicted if the mech anism of loss- 
cone r efilling is collisionless. To date, the iBerczik et al] 
(|2006( l simulations are the only ones that have success- 
fully followed the evolution of a binary SMBH all the 
way from galactic to sub-parsec scales, and that can ro- 
bustly be scaled to real galaxies due to the lack of an 
appreciable iV-dependence of the e yolution rates. 

The final binary separation in the IBerczik et al.l (j2006f ) 
simulations was ~ 0.05ah, just small enough that gravita- 
tional wave emission would induce coalescence in 10 Gyr 
(assuming a scaling that assigns a mass IO^Mq to the bi- 
nary SMBH). In this paper, we repeat the Berczik et al. 
simulations using a new iV-body code that incorporates 
the post-Newtonian {VN) corrections to the equations of 
motion of the SMBHs. This allows us to follow the evo- 
lution of the binary all the way to final coalescence. In 
addition, we are able to estimate the strength of the grav- 
itational waves emitted during the final inspiral. A key 
result is that the eccentricity of the binary can remain 
high until shortly before coalescence. Both the strength 
of the emitted gravitational waves, and the timescale for 
coalescence, are strongly dependent on the binary eccen- 
tricity. Highly eccentric black hole binaries would repre- 
sent appropriate candidates for forthcoming verification 
of gravitational radiation through the planned mission of 
LISA. 

The outline of this paper is as follows: in Section 2 we 
give a description of the numerical methods and initial 
models used in this work. In Section 3 we then present 
the results of a set of A^-body simulations using both 
the classical, i.e., Newtonian gravity (Sec. 3.1) and the 
relativistic, i.e., post-Newtonian gravity (Sec. 3.2). In 
Section 4 we discuss our results in a broader astrophysi- 
cal context, focusing on the SMBH merger as sources of 
gravitational wave emission. In Section 5 we finish with 
the conclusions. 

2. NUMERICAL MODELING 

For the A-body simulations presented in this work 
we use a mo dified version of th e publicly available'* (/3- 
Grape code (jHarfst et al.ll2007f). This cod e is based on 
an NBODYl-fike algorithm ijAarsethI 119990 . including a 



hierarchical timest ep scheme and a 4*^-orde r Hermite in- 
tegrator (see, e.g.. lMakino fc Aarsethl[T99^ . The gravi- 
tational acceleration a and its first time-derivative da/dt 
(also called: jerk) between the field particles - represent- 
ing the 'stellar' component in the galactic nucleus - are 
calculated using the special-purpose hardware Grape- 
6A (Fukushige, Makino & Kawai 2005). We also use the 
Grape hardware to calculate the interaction between 
the field particles and the black holes (and vice versa). 
The acceleration between the two black hole (BH) par- 
ticles, the corresponding jerk as well as the relativistic 
corrections, are all calculated using the hosts CPU. 

To determine the timestep bin fo r each individu al par- 
ticle we use the standard criterium (|Aarsethlll985[ ) of the 
form: 



,(2) I 



,(3) I 



,(2), 



(2) 
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where , a.; , and a,- are the acceleration, its first and 
fc'^ time-derivatives, respectively, of a given particle i. 
For details on how t o obtain the high er-order time- 
derivatives we refer to lHarfst et al.l (|20G7[ ). Whereas we 
use 77 = 77^ = 2 X 10^^ for the field particles, we use a 
smaller 77BH = 0.1?;^ = 2 x 10"'^ for the two BH par- 
ticles in order to guarantee an adequate conservation of 
energy and momentum (in the Newtonian gravity sim- 
ulations). Moreover, the two BHs are always advanced 
synchronously using the smaller AisH of the two. 

The simulations have been carried out on the dedi- 
cated high-performance GRAPE-6A clusters at the As- 
tronomisches Rechen-Institut in Heidelberg,^ at the Ro- 
chester Institute of Tcchnoloy,^ and at the Main Astro- 
nomical Observatory in Kiev.^ 

2.1. Initial Conditions 

The initial conditions of the A^-body models used in 
this work are chosen to be similar to the one s used by 
Berczik et al. (|2005f ) and IBerczik et al.l (|2006[) . allowing 
for a direct comparison to their simulations and for an 
accordant interpretation of our results. Here, we briefiy 
describe the main aspects of the initial model set-up and 
refer the reader to the latter two publications (and ref- 
erences therein) for further details. 

The initial field particle distribution - representing 
the galactic nucleus - is set-up following the phase- 
space density distribution o f a King mode l ('Kina"l966r) 
with internal rota tion fe.g.. iLongaretti fc L agoutc 1996| 
lErnst et al.l l2007l and references therein). The mod- 
els have been set up with a central concentration pa- 
rameter of Wq ~ 6 and an initial rotation of ujq = 
y/9/ {4tt G Pc) f^o = 1-8, where G and pc are the gravi- 
tational constant and the central concentration, respec- 
tively. The net rotation in our models is chosen to be 
counter-clockwise, with the angular momentum vector 
being aligned with the z-axis. We adop t the standard 
A^-body units (see, e.g.. I AarsethI l2003bl and references 
therein) to our numerical models by setting both the 
gravitational constant G and the total mass A/* of the 

^ GRACE: see |http: / /www.ari.uni-heidelberg.de/gracej 
^ gravitySin iulator: see | ht tp:/ /www.cs.rit.edu~ grapecluster| 
golowood: |http://www.rnao^^i^v!ua7goTowoon7eng] 
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iV-body system to unity, and by setting its total energy 
to = -1/4. 

In this work the galactic nucleus is represented with 
total numbers of A^^, = 25 x 10'^ and 50 x 10^ field par- 
ticles, respectively, each with nine different random real- 
izations. We restrict ourselves to these relatively small 
particle numbers as a trade-off for the large set of simula- 
tions that we h ave p erformed in grand total. However, as 
iBerczik et al.l (|2006f ) have shown, the hardening rate of 
SMBH binaries in such rotating galaxy models is in this 
regime essentially independent of the field particle num- 
ber Ni,. The main difference between our models and 
the latter ones is the 7W treatment of the two SMBHs. 
Therefore the results of our simulations are expected to 
be quasi- TV-independent as well, since we use a rotation 
parameter of > 1-2 dBerczik et al.ii2006i) . The field 
particles have all equal mass rrii, = Mi,/Ni,. The gravi- 
tational softening length for the field particles has been 
set to = 10"'* (in model units). The same softening is 
also used for the star^BH interaction. The two BH par- 
ticles themselves are evolved without any gravitational 
softening, i.e., with esH = 0.0. 

We set the masses mi and m2 of the two BH particles 
to be one per cent of the nuclei mass M^,, i.e., in our 
case mi = m2 = 0.01. The two BH particles are ini- 
tially placed on the y-axis at 7/1^2 ~ ±0.3, respectively, 
within the z = plane. The BHs are given an initial 
velocity Vx, corresponding to roughly ten per cent of the 
local circular velocity, as derived from the enclosed mass 
of the underlying field particle distribution. With this 
we get roughly vi^2 ~ ±0.07, in our model units. Note 
that in this configuration the two black holes are initially 
unbound with respect to each other. 

To scale our A^-body results to real galaxies we consider 
the total energy of the system 



E 



GIVP 
R 



(3) 



where G, M are the gravitational constant and the total 
mass of our model system (typically some fraction of the 
bulge and cusp components in a galactic nucleus), and R, 
a are model dependent quantities, giving a scaling radius 
and a numerical constant of order unity, respectively. As 
an example, for a non-rotating King model with Wq = 6 
(which is used in our models as an initial configuration) 
we w o uld h ave R ~ rc as the core radius as defined by 
IkImI ([l96l and a = 0.0759. 

With G = 1 we can freely choose two of the three scales 
(energy, radius, or mass). In standard A^-body units the 
unit of mass is the total mass, i.e. M = 1, and the 
unit of energy is AE {E = —0.25). Therefore the radial 
unit is determined by the condition R = 4a (this value 
determines the value of R in A^-body units, e.g., for a 
Plummer model we have a = 37r/64 and R = 37r/16). 
The physical units of mass, length, energy, velocity and 
time are then given as: 



[L]- 
[E] 

[V] 



-M 
R 
4a 

= 4a ■ 



GAP 



R 



[T] = 



2V^ 
1 



GM 
R 

R' 



1/2 



1/2 



8a3/2 \GM^ 
The speed of light c in A^-body units is then 



(4) 



c = co/ [V] = Co • 4a 



GM 



R 



-1/2 



= 457 



M 



10" Mf^ 



'1/2 



/ R 



Vl03pC 



1/2 



(5) 



where cq is the speed of light in physical units, and for the 
second expression we have inserted a = —0.25. This il- 
lustrates that in the post-Newtonian regime the A^-body 
problem is not scale-free anymore - to fix the scale for c 
means fixing the radial scale or vice versa. In our sim- 
ulations we have chosen different values for the speed of 
light, i.e., c = 447, 141,44, and 14, in order to enhance 
the effect of the relativistic corrections. At the same 
time, this variation of c scales the Schwarzschild radius 
i?BH = 2GAfBH/c^ of the two BHs to values of 10""^, 
10^^, 10^^ and 10~'*, respectively, in model units. 

2.2. Relativistic treatment of compact binary systems 

Super-massive black hole binaries are expected to be 
subject to the emission of gravitational waves. Since the 
gravitational waves drain energy and angular momentum 
from the binary, its orbital elements will change in the 
course of their dy namical evolution. Alread y mo re than 
four d ecades ago iPeters fc Mathew"i ()1963f ) and iPetersI 
()1964f ) derived the change of energy E and of the orbital 
elements for a Keplerian orbit, namely its semi-major 
axis a and eccentricity e, under the influence of the gravi- 
tational wave (quadrupole) emission. The orbit averaged 
expressions for two compact objects with masses mi and 
m2 orbiting each other are given as: 



/da 
\dt 



de 
'dt 
dE 
'dt 



64 G^ mi 77i2(mi -I- m2) 



304 G"* mi m2(mi 



b a'' 
^3 , 



fie), 



- m2) 



15 c5a4(i„e2)5/2 
32 G"^ mi m|(mi + m2) 



121 
304* 



/(e), 



(6) 
,(7) 
(8) 



where the so-called enhancement factor /(e), given by 
the following function, shows a strong dependence on the 
orbital eccentricity e: 



/(e) 



1 



-7/2 



1 



73 
— f 

24 



37 

96*^ 



(9) 



We note that the rate at which energy E is lost due to 
the emission of gravitational waves depends strongly on 
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the eccentricity e of the orbit. Hence, compact binaries 
with high eccentricities (i.e., small a and e 1) are ex- 
pected to be strong sources of gravitational waves. The 
equations of Peters & Mathews describe the shrinking 
of the semi-major axis (corresponding to an inspiral of 
the two ob jects) and the circul arization of the orbit. Ac- 
cording to IPeters fc Mathews! (|1963f) and iPetersI (|1964f ) 
the typical timescale of coalescence due to the emission 
of gravitational radiation is given by 



64 G^mi TO2(mi -I- 7712) /(e) 



(10) 



where Cgr denotes the characteristic separation for grav- 
itational wave emission. 

To account for the relativistic effects in our numeri- 
cal A'^-body simulations we use the post-Newtonian for- 
malism. The equations of motion for a comp act bi- 
nary system are written in harmonic coordinates (jSchut j 
[19851) and are defined in the inertial A^-body frame of 
reference. As they are relativistic, they are (1) invari- 
ant under T^A/^expandcd Lorentz transformations, (2) re- 
duce to the geodesies of the TW-expanded Schwarzschild 
metric in the limit in which one of the masses goes 
to zero and neglecting the spin, and (3) are conser- 
vative if the 2. SPA/" radiation reaction term is turned- 
off. In fact, an isolated 27W' binary should conserve 
its generalized VAf integrals of m otion up to 0(l/c^) 
(jAndrade. Blanchet fc Favd [200lf ) . We implement the 
VAf equations of motion formulated in t he absolute Eu- 
clidean space and time of Newton (e.g., iBlanchetl l2006l . 
equation 168 therein). In the present work we apply all 
T^A/" corrections up to the order 0(l/c^), i.e., the 2.57W 
corrections is the highest order that we take into account. 
While the 2. 57W correction accounts for the emission of 
gravitational waves, the IVAf and 27W terms conserve 
energy. However, the latter two terms result in preces- 
sion of the orbital pericenter. Similar to the equations of 
motion in the center of mass frame (Kupi, Amaro-Seoane 
& Spurzem 2006), one can write the accelerations, e.g., 
for particle 1 of the binary, in the following form: 



ai 



Gm2 



[(1 +^)ni2 +6V12] 



(11) 



where ri2 is the separation of the two particles, ni2 the 
normalized relative position vector, and V12 is the rel- 
ative velocity. The two functions A and B contain the 
different orders of the 7W expansion and can be written 



■ Bi-pjsf 



■ B2VN' 



■ B2.5VM + O 



•(13) 



In this notation, the first post-Newtonian correction 
(17W) is, for example, given as: 



5Gmi AGm2 3, ,n , 

+ + :^(ni2-V2)' - ^rl- 

ri2 ri2 2 



+ 4(vrV2)-2v2] , 
SimA = 4 (ni2-vi) - 3 (ni2-V2) . 



(14) 
(15) 



We forgo from writing down the higher order 7W correc- 
tions due to their lengthy form - especially the one of 
the 2y W" correction. Th e co mplete expre s sions are given 
in, e.g.. IBlanchetl (|2006f ) and lKupi et all (|2006f ). 

Our numerical implementation of the VM corrections 
to the accelerations - and particularly to the correspond- 
ing jerk as required for the Hermite integration scheme - 
has been thoroughly test ed against o t her av ailab le codes, 
such a s the ones used bv lKupi et al.l (|2006( l and I AarsethI 
(|2007f ). A detailed descripti on of the numerica l tests and 
the results are given in (B ercntzen et al.ll2008[ ). An alter- 
native implementation specially tailored to the treatment 
of single massive BH within galact ic nuclei can be found 
in lLockmann fc BaumgardB ()2008t ) . Our implementation 
of the TW corrections allows us to turn on and off the dif- 
ferent TW orders separately at will. In the current set of 
models presented here we apply the selected 7W terms to 
the two BH particles at all times during the simulations. 

Within the 7W framework, the energy l ost due to grav- 
itational wave emission is given by (see IBlanchetl [20061 
equation 171 therein): 



dE _A G^mlm2 
dt 5 c^r\2 



(V1V12) -vfa + 2 



Gmi 



, Gm-2 



+ (ni2-vi)(ni2-vi2) 

(1^2)+0(c-^). 



3vi2 - 6 



ri2 

Gmi 
+ 

ri2 



-8- 
ri2 

52 Gto2 

y ri2 



(16) 



Note that this latter equation gives - in contrast to the 
orbit averaged expression in Eq. [5] - the instantaneous 
energy loss due to the emission of gravitational waves. 
We will apply both equations, i.e., Eq. [8]and Eq. [161 for 
comparison in our analysis. 

3. RESULTS 

3.1. Purely Newtonian simulations - the fiducial case 

In this section wc describe the results of our purely 
Newtonian simulations, i.e., classical models without any 
corrections. For these fiducial models, the evolution 
of the binary BHs and the field particles i s very similar 
to the corresponding models reported by iBerczik et al.l 
(j2006[ ). Due to the relatively high degree of rotation in 
the models used in this work, the initially axisymmetric 
nucleus models become dynamically unstable and rapidly 
develop a rotating triaxial (bar-like) structure. As a re- 
sult of this global instability and due to the dynamical 
friction against the stellar background, the two BHs are 
quickly funneled towards the density center of the nu- 
cleus, where they form a binary system in our models 
typically after some At sa 20, or some 30Myr. 

It actually turns out that certain binaries character- 
istics, such as the time when the binary forms (tform), 
as well as its orbital elements, are very sensitive to the 
initial conditions of the stellar distribution. Since the ec- 
centricity of the binary is a crucial quantity for its sub- 
sequent evolution towards relativistic inspiral and coa- 
lescence, it is important to obtain a statistical sample. 
Therefore we decided to use nine different realizations 
of the same galaxy nucleus model, by sampling the un- 
derlying distribution function with different initial ran- 
dom seeds. This way, e.g., we find eccentricities of the 
binaries typically in the range between 0.4 up to 0.99. 
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Fig. 1. — Example of the typieal evolution of the SMBH binaries semi-major axis a (left panel), its eccentricity e (middle panel) and its 
total energy E (right panel) as a function of time in the classical (Newtonian) gravity models. The lines show the average eccentricity after 
the binary has formed and the linear fit to the energy loss due to super-elastic scattering. All quantities are given in model units. 
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Fig. 2. — Average rate of energy loss due to purely Newtonian, 
stellar-dynamical effects as a function of the binaries eccentricity. 
The symbols represent models with 25k (open circles) and 50k 
(filled circles) field particles, respectively. The horizontal lines in- 
dicate the mean values of \dE/dt\ for our set of 25fc (full line) and 
50fc (dashed line) models. The mean value of the combined set is 
plotted with a dotted line. Note that in our models the energy loss, 
or equivalently the binaries hardening rate, is basically independent 
of the eccentricity and — in contrast to spherically symmetric nu- 
clei models - also essentially A^-independent within the indicated 
error-bars. 



The binaries in our simulations tend to form with the 
higher eccentricities around 0.9. However, it cannot be 
decided at this point whether stochastic encounters be- 
tween the BHs and the stars, or the global (non-linear) 
bar-instability plays the dominant role in determining 
the final binary properties. 

In Fig. [T] we show an example for the typical time evo- 
lution of the (Keplerian) orbital elements for one of the 
SMBH binaries in our 50k Newtonian models. In the 
panels, from left to right, we show the evolution of the 
semi-major axis a, the eccentricity e, and the energy E 
of the binary, respectively. Note that the extremely large 
values of a and e at early times arc due to the fact the 
the two SMBH are still unbound. Only when they have 
formed a gravitationally bound pair, i.e., when the bi- 
naries energy E becomes negative, the eccentricity e re- 
mains below a value of 1. As one can see in the middle 
panel, the eccentricity varies only mildly in the subse- 
quent evolution, i.e., when the binary is hard and in- 



FlG. 3. — Characteristic timescales as a function of the semi- 
major axis for a) the Newtonian hardening (thick full line) and 
b) Peters & Mathews (dotted lines) for orbits with eccentricities e 
ranging from 0.5 to 0.9 (from top to bottom) in steps of Ae = 0.1. 
The speed of light has been chosen as c ~ 500 in model units. 



teracts with the surrounding field particles mainly by 
super-elastic three-body scattering. Therefore, we calcu- 
late some mean eccentricity e for each model as indicated 
by the horizontal line in the middle panel of Fig. [TJ We 
find e to be a useful quantity to characterize the binary, 
but it is important to bear in mind that it provides only a 
first order approximation, since the real eccentricity may 
in fact vary with time due to stochastic encounters with 
field particles. 

We find that the loss of energy (Fig. [TJ right panel) due 
to super-elastic three-body encounters is almost constant 
after roughly t > 30. The straight line shows the result 
of a linear least-squares fit to E{t) after the binary has 
formed. In Fig.[2]we show the resulting slopes \dE/dt\ as 
derived from such linear fitting as a function of (1 — e^) 
for our different Newtonian models. * The horizontal lines 
in Fig. [2] indicate the mean values of dE/dt for both the 
separate and combined sets of our Newtonian 25k and 
50k models. 

Since E cx (1/a) for a Keplerian orbit, the rate dE/dt 

* This quantity (1 — e^) is used rather then e itself, because it 
enters with a power of —7/2 in the enhancement factor f(e) (see 
Eq.[9]in the text). 
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Fig. 4. — Time-derivatives of the orbital elements a, e and E (from left to right). The gray curves show the discrete differences as 
calculated from the model directly, and the black curves show the derivatives based on the Peters & Mathews formalism. The dotted curve 
in the right panel shows the results from the post-Newtonian expression (see Eg. I16l in the text). 



provides a direct measure for the hardening rate of the 
binary djdt (1/a). We find that the hardening rate shows 
an essentially A^-independence of both the eccentricity of 
the binaries' orbit and the number of field particles used 
in the models. 

The latter resu lt is in good agreement with 
iBerczik et al.l (l2006l ) who found only small variations 
in the hardening rate for models with particle numbers 
ranging from 25k to IM field particles, which is inter- 
preted as a result of an efficient loss-cone refilling as 
found in such triaxial, rotating nuclei models. 

For our two Newtonian sets of 25k and 50k models, we 
find mean values for \dEldt\ of about (1.01±0.19) x IQ-^ 
and (0.95±0.17) x 10"^, respectively. The hardening rate 
in our models is found to be essentially iV-independent 
within the given error margi ns. This find i ng is clearly 
in accord with the results of 'Berczi k et al.l ([2006'). The 
mean value of dE/dt for the combined set of models thus 
results in being roughly E = dE/dt w —0.001. 

One can use the results of the previous paragraph to 
make some simple estimates regarding the binary evolu- 
tion timescale in the purely Newtonian regime described 
above. For a Keplerian orbit we get: 



da _ 2a^ dE 
dt Gmi 7712 dt 



(17) 



We define the binary hardening time in the usual way as 



Thard = 



1 da 


-1 


1 dE 


a dt 




E~dt 



G 



mi 7712 



2a 



dE 



dt 



(18) 



Writing the rate of change of the binary energy in A^- 
body units as \dE/dt\ = 10""^ x A', where if sa 1, this 
becomes in physical units 



Thard = 9.80 X WK 



5/2 



G^inili^a 



A ctn in4 mi 7712 

4.62 X 10 yr- 



(lO8M0)'a VlOV/ V10"M 



5/2 



-5/2 



This expression shows clearly that a constant rate of en- 
ergy loss corresponds to a gradually increasing timescale 



for binary hardening. If we assume that the expression 
holds for arbitrarily small a, we can predict a time to 
coalescence of 



icoal — tf) ~ ■''hard(ao) 



- 1 



(20) 



^coal 

where ao — a(to) and acoai = a(tcoai) (jFerrarese fc FordI 
I2005f) . Setting ao = 1 pc and acoai = (10"^ . . . 10"'') ao 
gives coalescence times ranging from ~ 10 Myr to 1 Gyr 
for a galaxy with Mgai = IO^^Mq and = 100 pc. Thus 
we see clearly that the binaries in our galaxy models 
would have no difficulty reaching very small separations 
in a Hubble time, even in the absence of relativistic en- 
ergy loss. 

In Fig. [3] we plot both the calculated timescale Thard 
for the stellar dynamical hardening, as well as the rela- 
tivistic coalescence timescale tg^ (see Eq. [T0| for orbits 
of different eccentricity as a function of their semi-major 
axis. As one can see, the two involved timescales are 
of the order of only some hundred time units or Myr, 
respectively, for the range of eccentricities found in our 
simulations. These timescales are short enough to allow, 
on average, for the SMBH binaries to coalesce by the 
time the next galaxy merger would occur. 

In order to test whether the purely Newtonian stellar 
dynamics in our models is sufficient enough to bring the 
SMBH binary into the regime where relativistic effects 
start to become important, we have done the following 
test: we first calculate the discrete time-derivatives, i.e., 
finite differences, of the orbital elements and the total 
energy as directly computed from our simulations. An 
example is shown in Fig. 3] (gray curves). In a second step 
we can then use Peters & Matthews' formalism (Eqs. 2- 
4) in order to calculate the corresponding orbit-averaged 
rates, which encode the binary's secular drift due to the 
emission of GWs at the lowest, quadrupolar, order. The 
resulting curves are also plotted in Fig. |4] (black curve). 
Finally, we use the post-Newtonian expression for the (in- 
stantaneous) energy loss due to gravitational wave emis- 
sion (Eq.[l6l). 

One can see that by the time t w 100 in model units 
(or some 150 Myr), the two black holes have reached a 
(Sep^faltion that is small enough for the relativistic ef- 
1 fe?^^^ ^^'^ reach the same order of magnitude as 

the (Newtonian) perturbations from the field particles. 
This analysis demonstrates that with these rotating King 
models, the SMBH binary quickly (within a cosmological 
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Fig. 5. — Time evolution of a SMBH binary orbit in a Newtonian + 2.5'PA'^model (using c = 447). Layout for the upper and lower panels 
are the same as in Figs. [T] and U] respectively. 
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Fig. 6. — Example of the time evolution a binary in the Newtonian + full PA^ models with c = 447. Same layout as in Fig. [5] 



context) hardens to a point where it reaches the relativis- 
tic regime and the final parsec problem can be overcome. 
In the next section, we study the self-consistent evolu- 
tion of the (non-spinning) SMBH binaries in the field of 
stars of the galactic nuclei remnant, including the post- 
Newtonian corrections up to 2.57Worder. 

3.2. Models including Post-Newtonian Corrections 

To test if the final parsec problem can really be over- 
come sclf-consistently just by stellar-dynamical effects 



and if the binary black holes in our models eventually 
reach the relativistic inspiral regime, we take the follow- 
ing approach: 

We start with models having exactly the same initial 
conditions for the set of 25k and 50k models described in 
the previous section, but this time we take the 7W correc- 
tions for the two black hole particles into account. Most 
numerical simulations of similar kind which are known 
to us from the literature are restricted ~ for simplicity 
- to the first dissipative TWterm in the expansion, i.e., 
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the 2. Sm/" correction. The effects of the IPTVand 2VM 
, which both are conservative and are known to result 
in a pericenter shift, have often been neglected. How- 
ever, these corrections are of lower order in (w/c) as com- 
pared to the 2.57Wand thus should clearly affect the bi- 
nary evolution. Therefore, we decided to run all sets of 
simulations presented below with (a) only including the 
2. 57W correction and (b) including the full PA/" correc- 
tions up to 2.hVM. Since higher order corrections such 
as the STWtcrm in the equations of motion arc of order 
c^^ or smaller, their effect on the dynamical evolution of 
the BHs are expected to be small as well, and will not 
affect the main results presented in this work. Hence, 
we will loosely refer to the models in our set (b) as the 
ones with jull 7W corrections. We should note here that 
the TW corrections in our runs are applied permanently, 
independent of the BH s eparation or t heir relative veloc- 
ities (see for comparison I Aarsethll2007[ ) . Finally, in order 
to enhance the relativistic effects we use different values 
of the speed of light c (in models units). We note that by 
fixing the values for G and c, the ratio between the units 
of mass and length is also fixed. Therefore, one has only 
one free parameter for setting the corresponding physical 
units in contrast to classical A^-body simulations, which 
are scale-free in that respect. 

In Figs. [5] and [6] (upper panels) wc show some examples 
of the evolution of the semi-major axis a, the eccentricity 
e and the total energy E of the binary, respectively. ^ The 
corresponding time derivatives, i.e., the discrete, Peters 
& Mathews and VNas defined in the previous section are 
shown in the lower panels of these figures, respectively. 
Both simulations shown in these figures start with iden- 
tical initial conditions and only differ by the orders of 
PA/" corrections taken into account. 

At early times, the dynamical evolution of the binary 
is qualitatively very similar as compared to the evolu- 
tion in our purely Newtonian simulations (sec Figs. [1] 
and HI): the two BH particles are rapidly driven towards 
the nuclei center by dynamical friction and the triax- 
ial bar-structure in the stellar nucleus. The black holes 
form a pair at roughly t = 2b (or some 37Myr, corre- 
spondingly) . Follow ing the results and interpretation of 
iBerczik et all (j2006f ). the initially full loss -cone around 
the black hole binary is quickly depleted and the evo- 
lution of the hard binary thereafter takes place in the 
empty loss-conc regime. Due to the internal rotation and 
the triaxial structure of the nucleus the loss-cone is re- 
populated constantly by centrophilic orbits with a rate 
which is found to be independent of the number of field 
particles in the model. 

During this phase the binary evolution is mainly dom- 
inated by just the Newtonian stellar-dynamical effects. 
The average rate of energy loss during this phase is again 
of the order of —0.001 as has been also found in the New- 
tonian simulations (Fig. [2|). Furthermore, this value is 
independent of the eccentricity with which the binary 
has formed and remains constant during the Newtonian 

® Note that, strictly speaking, the exact expressions for the 
semi-major axis a, eccentricity e and energy E would require VJ^ 
(i.e., VPM and 27W) corrections by their own (Memmesheimer, 
Gopakumar & Schafer 2004). However, for a direct comparison 
with the Newtonian models presented in the previous section, we 
stay with the classical expressions throughout the rest of this work, 
unless stated otherwise. 




50 100 150 200 

time 



Fig. 7. — Comparison of models using 2.57-iV corrections only 
(green lines) and using the full 'PA^ corrections (red lines). Shown 
are the inverse of the semi-major axis (thick lines) and the eccen- 
tricity (thin lines) as a function of time. 



dominated phase. That relativistic effects are not yet im- 
portant during this phase of evolution is supported in the 
lower panels in Figs. [5] and [6] in which wc compare the 
results from the simulations with the ones as expected 
by applying Peters & Mathews equations. 

The relativistic changes of the orbital elements become 
of the same order of the Newtonian ones after roughly 
i > 120 (in Fig.O and t > 60 (in Fig.©. More generally, 
we find that the time of transition from the Newtonian 
dominated regime to the one in which relativistic effects 
are of about the same order, depends strongly on the 
(mean) eccentricity e of the binary's orbit. This result 
is not surprising since as we have already stated before, 
the dissipation rate due to gravitational wave emission 
becomes already important at earlier times for more ec- 
centric binaries. 

Finally, at about times t > 140 and 70 for Figs. [5] and [6l 
respectively, the rate of energy loss increases significantly 
due to the emission of gravitational waves and becomes 
very non-linear. The loss of orbital energy and angular 
momentum due to the 2.57W corrections leads to the 
inspiral and circularization of the orbit, and the binary's 
dynamics in this phase is dominated by the relativistic 
effects, i.e., the binary almost fully completely decouples 
from the surrounding nucleus. 

It is noteworthy that these are the first systematic as- 
trophysical TV-body simulations which actually follow the 
evolution of the supermassive black holes from their un- 
bound state, i.e. with order kpc separation, down to the 
sub-parsec scale close to relativistic coalescence and thus 
overcome the final parsec problem sclf-consistently. It 
seems that the origin of the latter problem partly has 
been the use of over-simplified models for the galactic 
nuclei, i.e., assuming spherically symmetric distributions 
without (net) angular momentum. 

4. DISCUSSION 

We have presented one of the first extensive sets 
of stellar-dynamical A'^-body simulations including full 
post-Newtonian corrections up to 2.57W for the dom- 
inant two-body interaction between two supermassive 
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black holes, which covers self-consistently all evolu- 
tionary phases from an initially unbound SMBH bi- 
nary with separations of order 10'^ pc down to the fi- 
nal relativistic coalescence (but see Aarseth 2003a for 
an early pioneering work). Our initial conditions are 
still special but plausible: a "young" galactic merger 
remnant is represented by a flattened, rotating stel- 
lar nucleus with the two SMBHs at a separation of 
initially some 10'^ pc. This is a more general class 
of initial m odels than used in most other papers on 
the subject (iMakino et al.lll993tlMilosavlievic fc Merritt 



2Q01; Hemsexi dorf. Sigurdsson fc Spurzemll2002l: 



Aarseth 



2003ai iMakino fc Funatd l2004l : iBerczik et all |2005( ). 
This could be seen as one possible solution of the long- 
standing final parsec problem for black hole mergers 
(|Begelman et al.lll980D . as in this class of models we find 
convergence in the hardening rates at low TV. This is 
the result of a dynamical regime where the supply of the 
stars into the loss-cone region is essentially collisionless 
and presumably guaranteed by a family of centrophilic 
orbits associated with this non-spherical potential (cf. 
Sections 1.1 and 5). 

In the following we discuss in further depth some dy- 
namical properties of our models and observational con- 
sequences of our results for gravitational wave instru- 
ments. 

4.1. Effect of the different VAf orders 

We have performed of order 150 different simulations, 
varying the initial data (statistical realization, particle 
numbers) and different physical scalings (speed of light, 
different 7W orders). In all models that include the VM 
equations of motion for the black holes we are able to 
follow the binaries' evolution down to the relativistic in- 
spiral and virtually to coalescence. Based on our simula- 
tions we find that inclusion of the ITWand 27W correc- 
tions, which have been neglected in most ear li er wor ks 
(again with the notable exception of lAarsethl ()2003af )). 
strongly affects the dynamical evolution of the binary. 
Fig. [7] shows this, where we directly compare the evolu- 
tion of 1 /a and eccentricity as a function of time for two 
otherwise equal models with full TWand 2. 57W correc- 
tions only, respectively. The time of coalescence differs 
by almost a factor of two between these two cases. We 
find that part of the reason for this diverging behavior 
is the different eccentricity with which the binary forms. 
Using full TWcorrcctions rather than only 2.57Wchanges 
the details of the binary's trajectory already during their 
first encounters. Consequently later on, the binary forms 
with different orbital elements. Since all close encounters 
in the system - up to the one leading to the binary for- 
mation - are stochastic (due to the dynamical instability 
of close few-body encounters) we do do not find any pre- 
ferred systematic trend in our models towards higher or 
lower eccentricities between the two types of simulations 
as a function of the type of VAf corrections used. 

4.2. Orbital properties of SMBH binaries 

In Fig. [5] we show the distribution of eccentrici- 
ties e obtained from the sample of all our simula- 
tions. It clearly shows that the binaries form pref- 
erentially with very high eccentricities. This is fa- 
vored by th e transient triaxial feature in our mod- 
els (see also, IBerczik et al.l [20061) which brings the two 
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Fig. 8. — Distribution of eccentricities of the binaries in the New- 
tonian regime. The white histogram shows the normaUzed distri- 
bution of eccentricities of SMBHs in the Newtonian evolutionary 
phase. The gray histograms shows the distribution of eccentricities 
at later times, i.e., when the binary have reached a separation of 
100 Schwarzschild radii. Note that at these short separations the 
majority of the orbits still has not yet fully circularized. 



black holes together with relatively small impact pa- 
rameter. Since galactic mergers lead dominantly to 
the formation of ste llar bars or other triaxial config- 
urations (see, e.g. , iNa ab. Khochfar. fc Burkcrt' '2006!; 
I Jesseit et all 120071 : iJohansson. Naab fc Burkcrt 200a). 
we expect that high eccentricities in SMBH binaries after 
galactic mergers are a robust phenomenon occurring in 
many real cases. Furthermore, the merger of spherical 
models of galactic nuclei along nearly-parabolic trajecto- 
ries also leads to the formation of SMBH binaries with 
a similar d istribution of high eccentricities (jPreto et al.J 
Isubmittedl ). 

It is known that during the relativistic inspiral, cmis- 
sion of GWs leads to the circularization of the binary 
(|Petersl 11964 ). In the same Fig. [8l we therefore also 
show the (shaded) histogram for the binary's eccentric- 
ities when their semi-major axis has fallen down to 100 
Schwarzschild radii i?BH, for the runs with the realis- 
tic values of c = 447. We can clearly see that, although 
circularization has substantially reduced the binaries' ec- 
centricity, the remaining distribution at such small sep- 
aration is still significantly different from zero. Since 
the orbital eccentricity of the binaries is very important 
to predict gravitational waveforms from these binaries, 
and the expected distribution of their orbital parameters 
strongly determines the complexity of the data analysis 
for gravitational wa ve instruments suc h as the planned 
LISA satellites fe.g.. lBabak et al.ll2008[ ). we are currently 
generalizing the present study to a wider range of more 
realistic initial conditions.lt may be interesting to note 
that this type of statistics provides an astrophysically 
motivated set of initial conditions for numerical relativ- 
ity simulations of SMBH binary mergers. During the 
last couple of years several groups have made significant 
progress in modeling black hole mergers by the solution 
of the full Einst ein equations and of numerical-relativity 
simulations (cf. iRezzolla et al.l (|2008f ) and citations 8-13 
therein). 
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Fig. 9. — Example of the time a binary spends at different semi- 
major axis intervals. The solid line shows a Gaussian fit to the 
histogram used to determine expectation value of a and a. 
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Fig. 10. — Expectation value of (a) as a function of the orbital 
eccentricity for simulations with c = 447. The fitted line only takes 
the - more realistic - full T-I/V^ models into account. 



As we have described in Sec. 13.21 we can distinguish 
three dynamical phases of SMBH binary evolution: (1) 
Newtonian phase, (2) mixed phase and (3) relativistic 
phase. It is thus interesting to see how much time a bi- 
nary will spend the corresponding separations. In Fig. [9] 
we plot a histogram showing the time the binary spends 
at a given semi-major axis normalized to the total bi- 
nary lifetime, measured from its formation till coales- 
cence. The binary orbit of the model selected for Fig. [5] 
spends some ten per cent during its evolution at semi- 
major axis of 10""^, or 1 pc in physical units. We fit a 
Gaussian distribution to the binned data to provide a 
rough measure of the average value (a) and of the dis- 
persion around it. 

In Fig. [To] we then plot (a) as a function of the orbital 
eccentricity. While the SMBH binaries with high eccen- 
tricity tend to have larger (a), the fitted slope of the dou- 
ble logarithmic plot of (a) vs. (1 — e) is —0.87. We con- 
clude therefore that the average pericenter rp = (a) (1 — e) 
of the SMBH orbit at any given average (a) has a remark- 
ably small variation, i.e., is nearly constant. Therefore ~ 



Fig. 11. — Probability to find the SMBH binaries with separa- 
tions below some given maximum semi-major axis. The different 
lines correspond to orbits with different eccentricity as given in the 
legend. 



to first order - the orbits of SMBH binary in our models 
form a one-parameter family. This result is relevant and 
interesting for the determination of gravitational wave 
backgrounds in the ultra-low and low-frequency areas 
from SMBH binaries in the universe: Eccentric SMBH 
binaries emit gravitational waves with a spectrum of fre- 
quencies g{n,e), where n denotes the higher harmon- 
ics over the basic mode of the circular orbit (see Ap- 
pendix). The harmonic which leads to the maximal emis- 
sion of gravitational waves ninax(e) is for high eccentric- 
ities (e > 0.8) completely determined by the frequency 
of mo tion at pericenter nperi cx: ( 1 — e)""^/^ (jPierro et al.l 
120011: lAmaro-Seoane et al.|[2008[ ). Thus for every value 
of average semi-major axis (a) there exists a typical peri- 
center value and thus a typical gravitational wave spec- 
trum. The expected signals of these objects (in particular 
highly eccentric ones) lie in the pulsar timing and lower 
LISA band s (see for more details, incl i iding triple black 
holes , also lEnoki fc Nagashiinal l2007t iHofi'man fc Loebl 
[200l . 

Fig. [Tl] illustrates the cumulative probability to find 
the SMBH binary with a semi-major axis smaller than 
a. For each of our runs with full VN terms and real- 
istic value of c = 447 one curve is plotted in the fig- 
ure. From this information one can deduce that SMBH 
binaries with high eccentricity stay longer at larger a, 
which is an information consistent with the previous fig- 
ure. The probabilities are normalized for each run sepa- 
rately, therefore these data cannot be used to determine 
probability distributions of a for samples of SMBH bina- 
ries in the universe. 

4.3. Timescales for SMBH coalescence 

In Sec. 13.11 we estimated the time required for the 
coalescence of SMBH binaries in our models driven by 
stellar-dynamical effects only. On the other hand, using 
the pure relativistic estimate of the coalescence timescale 
tgr (Eq. [T0|) . we find for typical semi-major axis val- 
ues of ah ~ 10~^ . . . 10^^ and eccentricities of e « 
0.5... 0.99 that the resulting timescale covers a range 
of tgr w 10^ . . . lO^"' in our model units. The interpreta- 
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Fig. 12. — Merging time Tmerge of the SMBH binaries for models with differ ent values of c (c = 14:black; c = 44:red; c = 141:green; 
c = 447: blue. The shaded region indicate the predicted merging times using Eg. I21l as given in the text. The upper and lower boundaries 
of the shaded regions account for different formation times Tform of the binary. Since the Newtonian dominated regime lasts longer for 
models with c = 447, the eccentricity of the orbit varies stronger and results in a larger scatter in the plot. The are no clear differences 
between models with 2.57-iVonly and full 7W corrections. 



tion of these numbers is that for small semi-major axes a 
and high eccentricities the two SMBHs basically directly 
plunge, while for larger a and moderate e the relativis- 
tic timescale becomes unreasonable large. The evolution 
and hardening of the binary in such case would there- 
fore mainly be driven by the Newtonian stellar-dynamics 
over some time period. Eventually, the semi-major axis 
may reach small enough values for the emission of grav- 
itational waves to become efhcient. The binary harden- 
ing would then equally be driven by both the classical 
Newtonian and the relativistic effects, before the latter 
eventually takes over, i.e., when the binary dynamically 
decouples from the stellar background. 

To get a more precise estimate of the time until rela- 
tivistic coalescence, we numerically integrate the follow- 
ing differential equation: 

T_ = T,_+/[(|) +(|) 1.^,(21) 

J / hard V '^'^ / P&mJ 

where the corresponding da/ dt are given by Eq. [5] and 
1171 respectively. For the integration we start with semi- 
major axis separation of some 10~^ and different eccen- 



tricities. For the formation time of the binary we assume 
Tform = 15 and 35 as a lower and upper limit, respec- 
tively. 

In Fig. [T2] we plot the coalescence time as measured 
from our simulations for different values of c as a function 
of eccentricity. The results are then compared to the 
results obtained from integration of the above Eq. [^Ij 
Considering the different times it takes the binaries to 
form a bound pair in the simulations our results agree 
(within the error limits) with the theoretical value. 

Not surprisingly, we find that the merging time in- 
creases with increasing values of c, since the relativis- 
tic corrections decrease for simulations starting with the 
same configuration. However, since the hardening rate 
due to super-elastic scattering with the binary is found to 
be constant in our simulations, we expect an upper limit 
of Tmorg of roughly 500 time units. It is important to note 
that in all our simulations the two black hole coalesce in 
less then 0.5 Gyr after they have formed a bound pair! 
Therefore we conclude that our results provide a dy- 
namical substantiation to t he picture of pr o mpt S MBH 
coalescence advoca.ted b y iVolonteri et al.l ()2003D and 
ISesana et all (|2005ll2007t) . Moreover, we should also note 
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that these times for black hole coalescences are of the 
same order as the average interval between consecutiv e 
mergers (cf. e.g. Fig. 2 i n iKhochfar fc BurkertI (|200lh . 
and lKhochfar fc BurkertI (j2006f )). Such timescales make 
it unlikely that there is enough time for a third black hole 
to interact with the binary before it has coalesced. If 
three-body interactions between SMBHs were the norm, 
then it would be quite difficult to understand the notably 
small amount of scatter in the M — a relation, as well as 
the scarcity of observational evidence for off-centered nu- 
clei. Note, however, that some fraction of triple interac- 
tions of SMBHs may be present in the universe and could 
cause extremely large eccentricities (e.g. through Kozai 
oscillations) and become visible through pulsar timing or 
even in the lower LISA band at comparatively large sepa- 
rations (orbital time scales) ( Iwasawa. Funato &: Makinol 



Iwasawa. Funato fc Makinof 



2006t iHoffman fc Loebl [20071: lAmaro-Seoane et al.l l2008l: 



2008ir 



4.4. Gravitational Wave Signal 

Merger events of SMBH binaries are expected to be 
one of the brightest possible sources of gravitationa l 
wave em i ssion to be detected by LISA (|Danzmannlll997l : 
iPhinnevl |2005[ ): The Laser Interferometer Space An- 
tenna, a proposed space-borne gravitational wave obser- 
vatory, scheduled for launch in 2018+ , is most sensitive 
to gravitational waves in the low-frequency regime (10~^ 
- 0.1 Hz); for SMBHs of more than about lO^Mg LISA 
will pre ferentially detect t he final merger and ringdown 
phases (jBabak et al.ll200"8l ) , while for smaller masses ear- 
lier evolutionary phases of SMBH evolution will become 
detectable in the LISA band, in particular for eccentric 
SMBHs (see discussion and references in Sect. 14. 3p . 

In most of the present literature the strain amplitude 
of the gravitational wave is usually estimated under the 
assumption of circular orbits. This has been motivated 
by the idea that massive binary systems are expected to 
have been completely circulari zed due to the emission of 
gravi tational waves (see, e.g., lPeters|[r964l : ISesana et al.l 
|2005() by the time they enter the LISA frequency band. 
In this section we will use our results to test this assump- 
tion. The details on how to extract the information on 
the gravitational wave strain amplitude from our numer- 
ical simulations calculation are given in the Appendix. 

The space-based gravitational wave detector LISA will 
be sensitive to signals of inspiralling compact objects 
with up to the maximum total mass M12 ~ 7.5 x 10^ 
Mq, where M12 is the total mass of the binary. Due to 
the involved low frequencies of the signal, such sources 
are out of reach of current ground-based detectors. The 
physical units adopted in Sec. 12.11 translate into black 



hole masses of about 10^ M© in our simulations. Black 
hole binaries in this mass regime, however, will inspiral 
with frequencies even lower than those of the LISA band 
and thus would not be detectable by LISA before their 
actual coalescence. 

In order to follow the black hole binary inspiral up to 
frequencies in range of interest for LISA, i.e, / > 10"'^ 
Hz, we first need to rescale our models to binary masses 
of order of 10^ or 10^ Mq. As mentioned earlier, the 
speed of light in model units scales as c oc (i?/M)^/^, 
where R and M are the units for length-scale and mass, 
respectively. Since c is already given with a fixed value 
from our simulations, by changing our unit of mass 



by factors 10^^ — lO^'*, the length-scale automatically 
changes to 0. 1 or Ipc, respectiv ely. Finally, following 
the results of ISesana et al.l (|2007t ). we place the system 
at the redshift z = 4, where LISA's detection rate for 
objects of that mass range is expected to be significant. 

Before actually calculating the strain amplitude of the 
gravitational wave signal, we increase the time-resolution 
of our output data for the black hole's trajectory during 
the late phases of inspiral. This is done by evolving the 
binary in isolation, starting with initial conditions (in the 
binaries center of mass frame) as given from the large- 
scale simulations at a stage when it is already decou- 
pled from the surrounding stellar system, i.e., when the 
binary evolution is dominated by relativistic dynamics 
(compare Fig. [6]). We convince ourselves that the binary 
is really decoupled from the stellar system by compar- 
ing the evolution of its orbital elements to the one of the 
"low" time-resolution. In fact, during the late stages of 
inspiral the perturbations from stars vary over timescales 
much longer than both the orbital period and the time for 
coalescence (at that moment) which allows us to safely 
ignore them. 

Using the information of the combined time-series of 
the binary evolution we then cal culate the orbital ele- 
ment s in the VM generalization (jMemmesheimer et al.l 
l2004f l. These allow us to calculate the different mo des 
of the characteristic strain amplitude hc^n using Eq. IA6I 
from the Appendix. 

In Fig. [13] we show hen for the first six harmonics of 
the signal generated by the inspiral in one of our 50k 
models for (rescaled) binary masses of 2 x 10^ Mq (left 
panel) and 2 x 10'^ Mq (right panel). The LISA sensitiv- 
ity curve Sh{f) displayed in the figures was generated by 
the Online Sensitivity Curve Generator with default set- 
tings, corresp onding to a three year observation period 
(|Larsonl[2007l ). We confirm that the BH binaries in the 
mass range 10^ — 10^ Mq indeed enter the LISA frequency 
band during the relativistic inspiral in our simulations. 
The height of ft-c,n above the sensitivity curve in Fig. [13] 
provides an approximate indication of the signal-to-noise 
ratio (SNR; see also Appendix) and is found to be high 
enough for detection by LISA. 

As a consequence of the very high initial eccentricity 
with which the MBH binaries form in our simulations 
(compare Fig. [5]), they reach the LISA band with an ec- 
centricity which may be significantly different from zero 
- the exact distribution depending on the redshift z of 
the sources (Preto et al., in prep.). As a result, there 
will be a significant contribution to the measured power 
from higher harmonics /„ = n/orb/(l + z), with n > 2 
(note that ?7 = 2 for a circular orbit). This can be seen 
from the plots shown here, where several harmonics con- 
tribute significantly to the GW signal and are well above 
the threshold for detection. Furthermore, we find that, 
as the binary chirps in the LISA band, it circularizes 
quite rapidly with the consequence that the n = 2 har- 
monic always becomes dominant at the last stages of the 
inspiral. 

These findings may have an important impact on the 
accurate computation of the inspiral waveforms as well 
as potentially leading to an increased upper limit of the 
range of detectable masses by LISA. Based on our results 
we suggest that orbits with non-vanishing eccentricities 
should indeed seriously be considered for the LISA data 
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Fig. 13. — Locus of a SMBH binary from one of our simulations, for the first six harmonics, in the LISA sensitivity diagram during their 
final inspiral and coalescence. The black curve is the LISA sensitivity obtained from the online generator (see text). The dimensionless 
characteristic strain is plotted against the observed frequency. This case corresponds to a calculation including all PA'' corrections term 
up to the radiation reactio n at 2.57^ level. The orbita l parameters adopted in this case were the full 2'PM accurate tangential eccentricity 
et and semi-major axis a ijMemmesheimer et al.ir2004l '). Left panel: the binary MBH has total mass m\2 = 2 X 10^ Mq; right panel: 
mi2 = 2 X 1O®M0; placed at redshift z = 4. 



analysis. Further details and discussion of the conse- 
quences for the LISA detection of SMBH binary inspirals 
are beyond the scope of this work and will be published 
elsewhere (Preto et al., in prep.). 

5. CONCLUSIONS 

We present the first TV-body models which self-consis- 
tently follow the evolution of binary supermassive black 
holes in merged galaxies from kiloparsec separations 
down to gravitational-wave-induced coalescence, and in 
which the early evolution of the binary is driven by col- 
lisionlcss loss-cone rcpopulation, allowing the results to 
robustly be scaled to real galaxies. Our simulations in- 
clude post-Newtonian corrections to the equations of mo- 
tion of the SMBH binary up to order 2.^VN; wc show 
that inclusion of the energy-conserving WM and 2VN 
terms is also crucial for obtaining the correct time de- 
pendence of the binary orbital parameters. We identify 
evolutionary phases in which the binary is still evolving 
due to stellar encounters but in which relativistic correc- 
tions to its two-body motion are also important, thereby 
showing that the consideration of all TWorders in the A^- 
body simulations is necessary for an accurate prediction 
of its orbital elements evolution. The SMBH binaries 
in our simulations often form with large eccentricities, 
and these high eccentricities arc maintained during the 
Newtonian phases of the evolution. Wc show that the 
gravitational wave signal measured by an interferometer 
like LISA would contain significant contributions from 
high-order harmonics from GWs emitted from binaries 
at cosmological distances, and that the nature of the sig- 
nal can depend strongly on the dynamical history of the 
binary prior to its entering the gravitational wave domi- 
nated regime. 

We have shown that supermassive black hole binaries 
in galactic nuclei can overcome the stalling barrier and 
will reach the relativistic coalescence phase in a timescalc 
shorter than the age of the universe. A gravitational 
wave signal expected for the LISA satellite from these 
SMBH binaries is expected, in particular due to the high 



eccentricity of the SMBH binary when entering the rela- 
tivistic coalescence phase. We present for the first time a 
comprehensive set of models which cover self-consistcntly 
the transition from the Newtonian dynamics, dynamical 
friction phase (with yet unbound SMBH binary) to the 
situation when relativistic, post-Newtonian corrections 
start to influence the relative SMBH motion. After the 
shrinking time scale becomes very short the binary de- 
couples from the rest of the galactic nucleus and can 
be treated as a relativistic two-body problem. We follow 
this evolution formally to the coalescence of the two black 
holes using 7W terms of up to order 2.5PAAand deter- 
mine the gravitational wave emission in different modes 
relative to the LISA sensitivity curve. We find that the 
orbital parameters of a SMBH binary, when entering the 
LISA band, depends on the previous dynamical history 
- in particular there exists a phase where the SMBH bi- 
nary is still partially coupled to the stellar environment 
via three-body encounters, but relativistic corrections to 
its two-body motion already play a role. 

Since purely Newtoni an models of SMBH binaries in 
rotating galactic nuclei (jBerczik et al.|[2006l ) have shown 
a quasi-iV-independence of the stellar-dynamical driven 
hardening rates, we conclude that the results obtained in 
this paper regarding the time required until relativistic 
merger holds for galactic nuclei with realistic parame- 
ters. In this work we limited ourselves to models with 
particle numbers up to 50k only, however, due to the in- 
dependence of the results in the Newtonian phase from 
N in the current work and in Bcrczik et al. (200^, we 
do not expect any significant changes in our main results 
for models with N > 50k. 

It should be noted, however, that our results show a 
strong dependence on initial conditions, and that our 
initial model is a very simple approximation to a post- 
merger galactic nucleus. Therefore our conclusion from 
this work can only be that it is possible to reach relativis- 
tic coalescence in a reasonable time (10* years). Any pre- 
diction of event rates for LISA would, however, require 
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a more careful estimate of the distribution of parame- 
ters for a realistic set of mergers, e.g. b y taking data 
from s emi-analytic merger trees, similar to ISesana et al.l 
()2005f ). This is the subject of presently ongoing work. 

We have also examined the effect of using 2.57Walone 
or full corrections. Simple two-body TW experiments 
show significant differences of the BHs trajectories, and 
we find generally a dependence of the merging time in 
the two cases. Different orbital parameters of the SMBH 
binary in the final phase would have a major impact on 
the waveforms of the GWs. However, the computation 
of the latter also requires simulations with higher VAf 
corrections, at least up to order 0{c~^), i.e., 37Win the 
equations of motion. This is beyond the scope of the 
current paper. 

Results on whether our simulated SMBH binary will 
fall into the LISA frequency band and at which frequen- 
cies will be discussed in a cosmological context in a forth- 
coming paper. 
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APPENDIX 
APPENDIX MATERIAL 

Any gravitational wave carries energy: The total energy carried by a wave is E ^ ^if) where N{f) is the 
number of cycles the wave spends on a frequency interval A/ ~ / around the frequency /. It is customary to define a 
characteristic strain corresponding to an observation with duration r > N{f)/f by hc{f) = a/ N{f) h(f). In this 
case, the signal is not monochromatic and we observe its chirp as it shifts to higher frequency during the observation. 
However, in case the observation time r < N{f)/f, the signal is approximately monochromatic and its amplitude 
is essentially li mited by the o bserving time rather than by the intrinsic properties of the system and, as a result, 
/t. ^A/fffefn (lThornd ll987l). 

In the relativistic ins piral phase, the back-reaction from the GW emission (energy balance argument) dominates the 
binary's orbital decay ( Blanchetll2006f ). The number of cycles around a given frequency / can be estimated to be 

^(/) - /^//' - ^G-^/3A^-^/v^^/^ (Ai) 

where the orbit is described in terms of the (instantaneous) Kepler elements, the binary's chirp mass is M = 
'^^^r ru,^^ /(mi + 7712)^^^, and /r is the orbital frequency in the source's rest frame. Following iPeters fc Mathews! 



m 

(|1963[ ). the orbital frequency shifts at a rate given by 

The solution to Eq. IA2I blows up in a finite time Tcoai, which is used to denote the coalescence time. For a given 
frequency /r the time till coalescence can by calculated according to the following expression: 

-coal - ^ - ^G-/3A^-/3/-/3. (A3) 

Assuming a Keplerian orbital parametrization, the strain amplitude for the n'^ harmonic (after sky averaging) is 
given by 

37^2/3 £-5/3^5/3 



hnif) = ^hl + hi = — -^-^fV^,/^). (A4) 

Note that r{z) is the comoving distance (as a function of cosmological redshift z) to the source and g{n^ e) represents 
the normalized relative power spectrum in the n*^ harmonic of the GW signal (jPeters fc Mathew3ll963f ): 



9{n,e) = — 



2 

Jn^2{ne) - 2eJ„_i(ne) + - J„(ne) + 2eJ„+i(rie) - J„+2("e) 

n 



+ (1 - £.2) [Jn-2{ne) - 2J„(ne) + J„+2(ne)]' + ^^'("e) \ , (A5) 



where the J„ are Bessel functions of the first kind. The characteristic strain hc^n can therefore be written as 

A^</obsT 



07 r<5/S A//5/6 



(A6) 



37^2/3 q5/3j\^5/3 

hc,n if) = \/fr h ^J^^ Vain, e) ^/7 , N > f^bs t. 



Note that the frequency seen in the detector is given by /obs = /r/(l + z), where fobs is the frequency measured in 
the binary's reference frame. The knee frequency separates the two regimes 
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/knoc — ^ if) fobs 



1 / 5 cl5/8 



(1 + ^) 



(A7) 



Note that the knee frequency decreases if the observation time r increases. 

The GW signal is said to come into band when hc^nif) > (/ <S'/i(/))^^^, where Sh{f) is the instrumental strain 
noise spectrum (in Hz^^) and the brackets denote sky-averaging. The LISA sensitivity curve Sh{f) displayed in the 
figures was generat ed by the On line Sensitivity Curve Generator with default settings, corresponding to a three year 
observation period (jLarsonll2007[ ). The height of he above the sensitivity curve provides an approximate indication of 
the signal-to- noise ratio (SNR); however, the integrated SNR should be computed as 




(A8) 



